* ===============================

* Title: Replication files for 'Pitfalls in comparing Paris pledges'
* Project: Re-visiting Drummond et al. (2018)'s finding about the relationship between public opinion and GHG targets
* Author: Sam S. Rowan
* Contact: sam.rowan@nuffield.ox.ac.uk
* Date: 15 April 2019

* ===============================

* Organization

* 1. Load in targets and GHG data and reshape
* 2. Merge in Drummond et al. and Lee et al. data
* 3. Run main analyses
* 4. Make tables
* 5. Run auxillary analyses

* ===============================

**** Set working directory ****

** Set the working directory to the outermost folder containing 4 folders (code, data, figures, table)
cd "/Users/sr/Dropbox/Projects/DHSP/dhsp-july/"

* ===============================
**** Load in GHG and targets data; create variables for analysis ****
* ===============================
 
** Load in GHG data and reshape to wide

use "data/primap-1988panel.dta", clear
reshape wide primap_l_ghg primap_x_ghg, i(ccode) j(year)
order primap_x_ghg** primap_l_**, after(iso)

** Merge in targets data and reshape 
merge 1:1 ccode using "data/targets.dta", gen(_m1)

** Generate new "nominal" target variable, unconditional targets only: paris_headline
gen stdized_by = (100-target_paris_baseyear)/-100
gen stdized_bau = (100-target_paris_scenario)/-100
gen stdized_int = (100-target_paris_intensity)/-100

gen paris_nominal = stdized_int
replace paris_nominal = stdized_bau if stdized_bau != .
replace paris_nominal = stdized_by if stdized_by != .

drop stdized_by stdized_bau stdized_int

** Generate reference absolute emissions levels that emissions will be cut from under INDC targets
gen scenario_referencelevel = bau_emissionslevel*1000  /* put WRI data on the same scale as PRIMAP GHG data */
label var scenario_referencelevel "Reference absolute emissions level, scenario targets"

** GHG emissions for countries where LULUCF leads to negative values
* Find countries where ghg-lulucf is ever below zero because these will not fit the data well
egen min_ghg_l = rowmin(primap_l_ghg1988 primap_l_ghg1989 primap_l_ghg1990 primap_l_ghg1991 primap_l_ghg1992 primap_l_ghg1993 primap_l_ghg1994 primap_l_ghg1995 primap_l_ghg1996 primap_l_ghg1997 primap_l_ghg1998 primap_l_ghg1999 primap_l_ghg2000 primap_l_ghg2001 primap_l_ghg2002 primap_l_ghg2003 primap_l_ghg2004 primap_l_ghg2005 primap_l_ghg2006 primap_l_ghg2007 primap_l_ghg2008 primap_l_ghg2009 primap_l_ghg2010 primap_l_ghg2011 primap_l_ghg2012 primap_l_ghg2013 primap_l_ghg2014)
* list ccode country min_ghg_l if min_ghg_l <0
* Gabon, Romania, Bhutan, Tonga, Palau, Samoa, Fiji, Latvia, Montenegro, Chile, Bosnia and Herzegovina 

* Replace their values of lulucf, so that we use their ghg-ex values in later calculations
replace lulucf = 0 if min_ghg_l <0

* Which years to look for in coding reference levels?
table paris_referenceyear 

gen by_referencelevel = .  
* 1990 base years
replace by_referencelevel = primap_x_ghg1990 if paris_referenceyear == 1990 & lulucf == 0
replace by_referencelevel = primap_l_ghg1990 if paris_referenceyear == 1990 & lulucf == 1
replace by_referencelevel = primap_x_ghg1990 if paris_referenceyear == 1990 & lulucf == .
* 1994 base years
replace by_referencelevel = primap_x_ghg1994 if paris_referenceyear == 1994 & lulucf == 0
replace by_referencelevel = primap_l_ghg1994 if paris_referenceyear == 1994 & lulucf == 1
replace by_referencelevel = primap_x_ghg1994 if paris_referenceyear == 1994 & lulucf == .
* 2000 base years
replace by_referencelevel = primap_x_ghg2000 if paris_referenceyear == 2000 & lulucf == 0
replace by_referencelevel = primap_l_ghg2000 if paris_referenceyear == 2000 & lulucf == 1
replace by_referencelevel = primap_x_ghg2000 if paris_referenceyear == 2000 & lulucf == .
* 2005 base years
replace by_referencelevel = primap_x_ghg2005 if paris_referenceyear == 2005 & lulucf == 0
replace by_referencelevel = primap_l_ghg2005 if paris_referenceyear == 2005 & lulucf == 1
replace by_referencelevel = primap_x_ghg2005 if paris_referenceyear == 2005 & lulucf == .
* 2007 base years
replace by_referencelevel = primap_x_ghg2007 if paris_referenceyear == 2007 & lulucf == 0
replace by_referencelevel = primap_l_ghg2007 if paris_referenceyear == 2007 & lulucf == 1
replace by_referencelevel = primap_x_ghg2007 if paris_referenceyear == 2007 & lulucf == .
* 2008 base years
replace by_referencelevel = primap_x_ghg2008 if paris_referenceyear == 2008 & lulucf == 0
replace by_referencelevel = primap_l_ghg2008 if paris_referenceyear == 2008 & lulucf == 1
replace by_referencelevel = primap_x_ghg2008 if paris_referenceyear == 2008 & lulucf == .
* 2010 base years
replace by_referencelevel = primap_x_ghg2010 if paris_referenceyear == 2010 & lulucf == 0
replace by_referencelevel = primap_l_ghg2010 if paris_referenceyear == 2010 & lulucf == 1
replace by_referencelevel = primap_x_ghg2010 if paris_referenceyear == 2010 & lulucf == .
* 2012 base years
replace by_referencelevel = primap_x_ghg2012 if paris_referenceyear == 2012 & lulucf == 0
replace by_referencelevel = primap_l_ghg2012 if paris_referenceyear == 2012 & lulucf == 1
replace by_referencelevel = primap_x_ghg2012 if paris_referenceyear == 2012 & lulucf == .
* 2013 base years
replace by_referencelevel = primap_x_ghg2013 if paris_referenceyear == 2013 & lulucf == 0
replace by_referencelevel = primap_l_ghg2013 if paris_referenceyear == 2013 & lulucf == 1
replace by_referencelevel = primap_x_ghg2013 if paris_referenceyear == 2013 & lulucf == .
* 2014 base years
replace by_referencelevel = primap_x_ghg2014 if paris_referenceyear == 2014 & lulucf == 0
replace by_referencelevel = primap_l_ghg2014 if paris_referenceyear == 2014 & lulucf == 1
replace by_referencelevel = primap_x_ghg2014 if paris_referenceyear == 2014 & lulucf == .

replace by_referencelevel = . if paris_baseyear != 1 
label var by_referencelevel "Reference absolute emissions level, base year targets"

** Generate common INDC reference level
gen indc_referencelevel = scenario_referencelevel 
replace indc_referencelevel =  by_referencelevel if indc_referencelevel == .

** Generate target absolute emissions level under INDC
gen abs = 1-(abs(paris_nominal))

* For base year targets
gen by_targetlevel = by_referencelevel * abs
replace by_targetlevel = .  if paris_baseyear != 1
label var by_targetlevel "Absolute emissions target, base year targets"

* For scenario targets
gen scenario_targetlevel = .
replace  scenario_targetlevel = scenario_referencelevel * abs
replace  scenario_targetlevel = . if paris_scenario !=1
replace scenario_targetlevel = . if flag_target == 1 /* Check that no non-quantifiable targets get generated; use codebook flag_target to see labels */
replace scenario_targetlevel = . if country == "Mali" /* Mali has negative emissions target that the measures can't incorporate properly */
label var scenario_targetlevel "Absolute emissions target, scenario targets"

* Generate common INDC target level: merge these two variables
gen indc_targetlevel = by_targetlevel
replace indc_targetlevel = scenario_targetlevel if indc_targetlevel == .
replace indc_targetlevel = . if country == "Mali" /* Mali has negative emissions target that is poorly accounted for in the data */
label var indc_targetlevel "Absolute emissions target in INDC"

** Generate standardized emissions reductions, in terms of changes from different GHG years
* ===============================
forval yr = 1990/2014 {
	gen delta_`yr' = indc_targetlevel/primap_x_ghg`yr'
	replace delta_`yr' =  indc_targetlevel/primap_l_ghg`yr' if lulucf == 1
	replace delta_`yr' = (delta_`yr' -1)*-100 /* Direction: positive indicates emissions cuts */
	label var delta_`yr' "INDC target emissions as % of `yr' emissions"
}

** Emissions growth data for changes in GHG growth rate indicator

** Emissions for any particular year can be idiosyncratic/volatile so average them

* Early period average
egen avg_ghg_l_1992 = rowmean(primap_l_ghg1990 primap_l_ghg1991 primap_l_ghg1992 primap_l_ghg1993 primap_l_ghg1994)
egen avg_ghg_x_1992 = rowmean(primap_x_ghg1990 primap_x_ghg1991 primap_x_ghg1992 primap_x_ghg1993 primap_x_ghg1994)

* Late period average
egen avg_ghg_l_2012 = rowmean(primap_l_ghg2010 primap_l_ghg2011 primap_l_ghg2012 primap_l_ghg2013 primap_l_ghg2014)
egen avg_ghg_x_2012 = rowmean(primap_x_ghg2010 primap_x_ghg2011 primap_x_ghg2012 primap_x_ghg2013 primap_x_ghg2014)

* Historical emissions growth with the compound interest formula
gen ghg_growth_l_19922012 = (avg_ghg_l_2012/avg_ghg_l_1992)^(1/20)
gen ghg_growth_x_19922012 = (avg_ghg_x_2012/avg_ghg_x_1992)^(1/20) 

* Generate common historical emissions growth variable
gen hist_growthrate = ghg_growth_x_19922012
replace hist_growthrate = ghg_growth_l_19922012 if lulucf==1

** Generate growth rates for targeted emissions under INDC
gen indc_growth_l = (indc_targetlevel/avg_ghg_l_2012) ^ (1/18)
replace indc_growth_l = (indc_targetlevel/avg_ghg_l_2012) ^ (1/13) if paris_targetyear == 2025
gen indc_growth_x = (indc_targetlevel/avg_ghg_x_2012) ^ (1/18)
replace indc_growth_x = (indc_targetlevel/avg_ghg_x_2012) ^ (1/13) if paris_targetyear == 2025

** Generate common compliance growth rate variable 
gen indc_growthrate = indc_growth_x
replace indc_growthrate = indc_growth_l if lulucf == 1

* For BAU countries, we can see how their forecasted growth compares to their past growth
gen bau_growth_l = (scenario_referencelevel/avg_ghg_l_2012) ^ (1/18)
replace bau_growth_l = (scenario_referencelevel/avg_ghg_l_2012) ^ (1/13) if paris_targetyear == 2025
gen bau_growth_x = (scenario_referencelevel/avg_ghg_x_2012) ^ (1/18)
replace bau_growth_x = (scenario_referencelevel/avg_ghg_x_2012) ^ (1/13) if paris_targetyear == 2025

gen bau_growthrate = bau_growth_x
replace bau_growthrate = bau_growth_l if lulucf ==1 
replace bau_growthrate = . if target_paris_scenario == .

** Compare growth rates

* Targets versus historical
* ===============================
gen delta_indcgrowthrate = indc_growthrate - hist_growthrate 
replace delta_indcgrowthrate = -1* delta_indcgrowthrate /* Fixes the direction to match other indicators */

* Rescale paris_nominal to match Drummond et al. (ease interpretation)
replace paris_nominal = paris_nominal*-100 /* Fixes the direction to match other indicators */

* Rescale delta_growthrate because the values are really small (ease interpretation)
replace delta_indcgrowthrate = delta_indcgrowthrate*100


* ===============================
**** Merge in covariates data ****
* ===============================

** Merge in Drummond et al. 
merge 1:1 ccode using "data/Drummond, Hall, Sauer & Palmer CC 2018 Climate.dta", gen(_m2)
rename awareness drummond_awareness

** Merge in Lee et al. 
merge 1:1 ccode using "data/Lee et al. NCC 2015.dta", gen(_m3)
clonevar lee_awareness = aware
clonevar lee_serious = serious
drop wp5 region sub_region population2008 glo_tot efcon totbiocap va ps gove rq rl cc unaware rf_aware auc_aware top_aware top2_aware top3_aware ratio_aware not_serious rf_serious auc_serious top1_serious top2_serious top3_serious ratio_serious
* Countries are not in the data if they are not surveyed in Lee et al. 2015

** EU variables
summ co2emi gdpperus wgi spcias wpcias apcias lee_aware lee_serious if eu == 1

replace co2emi = 8.776641 if country == "EC"
replace gdpperus = 35167.71 if country == "EC"
replace wgi = .6101277 if country == "EC"
replace spcias = 241.1965 if country == "EC"
replace wpcias = .1532026 if country == "EC"
replace apcias = 51.63738 if country == "EC"
replace lee_awareness = 90.290001 if country == "EC"
replace lee_serious = 69.1277 if country == "EC"

* Rescale GDP into 1000s to ease interpretation in regression tables
gen gdp1000 = gdpperus/1000


* ===============================
**** Regression models ****
* ===============================

** Sample composition
reg unconditional co2emi gdpperus wgi spcias wpcias apcias lee_awareness 
gen model1 = e(sample)
reg paris_nominal co2emi gdpperus wgi spcias wpcias apcias lee_awareness 
gen model2 = e(sample)
reg delta_indcgrowthrate co2emi gdpperus wgi spcias wpcias apcias lee_awareness 
gen model3 = e(sample)
reg delta_2005 co2emi gdpperus wgi spcias wpcias apcias lee_awareness
gen model4 = e(sample)

** Common sample
gen common_sample123 = model3
replace common_sample123 = 0 if model2==0
replace common_sample123 = 0 if model1==0
tab common_sample123

gen common_sample23 = model3
replace common_sample23 = 0 if model2==0
tab common_sample23

save "data/dhsp-v4-july.dta", replace

** Replications
* ===============================

* Set working directory for tables
cd "tables"

** Table 1 ** 
* ===============================

** Model I: Original specification (Drummond et al. 2018)
* ===============================
bootstrap, reps(5000) seed(42): reg unconditional co2emi gdp1000 wgi spcias wpcias apcias lee_awareness
outreg2 using drummond_main, se tex dec(3) adjr2 10pct
bootstrap, reps(5000) seed(42): reg unconditional co2emi gdp1000 wgi spcias wpcias apcias lee_serious 
outreg2 using drummond_main, se tex dec(3) append  adjr2  10pct

** Model II: Re-coded headline variable
* ===============================
bootstrap, reps(5000) seed(42): reg paris_nominal co2emi gdp1000 wgi spcias wpcias apcias lee_awareness
outreg2 using drummond_main, se tex dec(3) append  adjr2  10pct
bootstrap, reps(5000) seed(42): reg paris_nominal co2emi gdp1000 wgi spcias wpcias apcias lee_serious 
outreg2 using drummond_main, se tex dec(3) append  adjr2  10pct

** Model III: INDC growth rate
* ===============================
bootstrap, reps(5000) seed(42): reg delta_indcgrowthrate co2emi gdp1000 wgi spcias wpcias apcias lee_awareness
outreg2 using drummond_main, se tex dec(3) append  adjr2  10pct
bootstrap, reps(5000) seed(42): reg delta_indcgrowthrate co2emi gdp1000 wgi spcias wpcias apcias lee_serious 
outreg2 using drummond_main, se tex dec(3) append  adjr2  10pct

* Set working directory for these outputs: get re-shaped in R and plotted
cd "../figures"

** Data for figure 3 ** 
* ===============================

** Coefficients for emissions reductions from different common years
* ===============================
* Bootstrap: each year takes about a minute
* Figures are built in R by reading in these tables, see figures code
forval yr = 1990(1)2014 {
	qui bootstrap, reps(5000) seed(42): reg delta_`yr' co2emi gdp1000 wgi spcias wpcias apcias lee_awareness 
	outreg2 using drummond_boot_awareness, keep(lee_awareness) stats(coef se ci_low pval ci_high) noaster dec(3) nopar append excel
	qui bootstrap, reps(5000) seed(42): reg delta_`yr' co2emi gdp1000 wgi spcias wpcias apcias lee_serious 
	outreg2 using drummond_boot_serious, keep(lee_serious) stats(coef se ci_low pval ci_high) noaster dec(3) nopar append excel
	dis `yr'
 }


* ===============================
**** Auxillary analyses and supplementary information ****
* ===============================

** Summary statistics
* ===============================
* Set working directory for tables
cd "../tables"

* Will need to combine these in Latex to make the summary statistics table in document
sutex unconditional if model1==1, labels minmax digits(1) file(summtables_unconditional.tex)
sutex paris_nominal if model2==1, labels minmax digits(1) file(summtables_nominal.tex)
sutex delta_indcgrowthrate if model3==1,  labels minmax digits(1) file(summtables_deltagrowth.tex)
sutex delta_1* delta_2* if model4==1, labels minmax digits(1) file(summtables_deltayears.tex)
sutex lee_awareness lee_serious co2emi gdp1000 wgi spcias wpcias apcias, labels minmax digits(1) file(summtables_covariates.tex)

** Regression models with common samples
* ===============================
* Set working directory for tables
cd "../tables"

* Model I: aware
bootstrap, reps(5000) seed(42): reg unconditional co2emi gdp1000 wgi spcias wpcias apcias lee_awareness if common_sample123==1
outreg2 using dhsp_commonsample123, keep(lee_awareness) stats(coef se ci_low pval ci_high) noaster dec(3) nopar append excel
bootstrap, reps(5000) seed(42): reg unconditional co2emi gdp1000 wgi spcias wpcias apcias lee_awareness if common_sample23==1
outreg2 using dhsp_commonsample23, keep(lee_awareness) stats(coef se ci_low pval ci_high) noaster dec(3) nopar append excel
* Model I: serious
bootstrap, reps(5000) seed(42): reg unconditional co2emi gdp1000 wgi spcias wpcias apcias lee_serious if common_sample123==1
outreg2 using dhsp_commonsample123, keep(lee_serious) stats(coef se ci_low pval ci_high) noaster dec(3) nopar append excel
bootstrap, reps(5000) seed(42): reg unconditional co2emi gdp1000 wgi spcias wpcias apcias lee_serious if common_sample23==1
outreg2 using dhsp_commonsample23, keep(lee_serious) stats(coef se ci_low pval ci_high) noaster dec(3) nopar append excel

* Model II: aware
bootstrap, reps(5000) seed(42): reg paris_nominal co2emi gdp1000 wgi spcias wpcias apcias lee_awareness if common_sample123==1 
outreg2 using dhsp_commonsample123, keep(lee_awareness) stats(coef se ci_low pval ci_high) noaster dec(3) nopar append excel
bootstrap, reps(5000) seed(42): reg paris_nominal co2emi gdp1000 wgi spcias wpcias apcias lee_awareness if common_sample23==1 
outreg2 using dhsp_commonsample23, keep(lee_awareness) stats(coef se ci_low pval ci_high) noaster dec(3) nopar append excel

* Model II: serious
bootstrap, reps(5000) seed(42): reg paris_nominal co2emi gdp1000 wgi spcias wpcias apcias lee_serious if common_sample123==1
outreg2 using dhsp_commonsample123, keep(lee_serious) stats(coef se ci_low pval ci_high) noaster dec(3) nopar append excel
bootstrap, reps(5000) seed(42): reg paris_nominal co2emi gdp1000 wgi spcias wpcias apcias lee_serious if common_sample23==1
outreg2 using dhsp_commonsample23, keep(lee_serious) stats(coef se ci_low pval ci_high) noaster dec(3) nopar append excel

* Model III: aware
bootstrap, reps(5000) seed(42): reg delta_indcgrowthrate co2emi gdp1000 wgi spcias wpcias apcias lee_awareness if common_sample123==1
outreg2 using dhsp_commonsample123, keep(lee_awareness) stats(coef se ci_low pval ci_high) noaster dec(3) nopar append excel
bootstrap, reps(5000) seed(42): reg delta_indcgrowthrate co2emi gdp1000 wgi spcias wpcias apcias lee_awareness if common_sample23==1
outreg2 using dhsp_commonsample23, keep(lee_awareness) stats(coef se ci_low pval ci_high) noaster dec(3) nopar append excel

* Model III: serious
bootstrap, reps(5000) seed(42): reg delta_indcgrowthrate co2emi gdp1000 wgi spcias wpcias apcias lee_serious if common_sample123==1 /*  */
outreg2 using dhsp_commonsample123, keep(lee_serious) stats(coef se ci_low pval ci_high) noaster dec(3) nopar append excel
bootstrap, reps(5000) seed(42): reg delta_indcgrowthrate co2emi gdp1000 wgi spcias wpcias apcias lee_serious if common_sample23==1 /*  */
outreg2 using dhsp_commonsample23, keep(lee_serious) stats(coef se ci_low pval ci_high) noaster dec(3) nopar append excel


** Tabulation of public awareness by target quantifiabilty 
* ===============================
gen observed = 1 if flag_target==0
replace observed = 0 if observed==.
ttest lee_aware, by(observed)
* Unobserved: 54.5 [49.7, 59.4]; Observed: 69.7 [64.1, 75.2]



* ===============================
* ===============================
* End of analyses
* ===============================
* ===============================




